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We study the Dirac equation for quasiparticles in gapped graphene with two oppositely charged 
impurities by using the technique of linear combination of atomic orbitals and variational Galerkin- 
Kantorovich method. We show that for sufficiently large charges of impurities the wave function of 
the occupied electron bound state of the highest energy changes its localization from the negatively 
charged impurity to the positively charged one as the distance between the impurities increases. 

This migration of the electron wave function of supercritical electric dipole is a generalization of the 
familiar phenomenon of the atomic collapse of single charged impurity to the case where electron-hole 
pairs are spontaneously created from vacuum in bound states with charge impurities thus partially 
screening them. 

PACS numbers: 81.05.ue, 73.22.Pr 


I. INTRODUCTION 

Atomic collapse or supercritical Coulomb center instability is a well-known phenomenon in quantum electrodynam¬ 
ics. For the electron in the Coulomb field of a nucleus of charge Ze, atomic collapse takes place for Z > 170 
when the lowest energy electron bound state dives into the lower continuum. This leads to the spontaneous creation 
of electron-positron pairs with the electrons screening the positively charged nucleus and the positrons emitted to 
infinity. Since supercritically charged nuclei are not encountered in nature, this phenomenon was never observed in 
QED. The situation essentially changes in graphene whose effective coupling constant e’^Khvp) ~ 2.2 {vp ~ c/300 is 
the Fermi velocity) exceeds unity that drastically decreases the value of the critical charge in graphene [2-13 . 

Experimental observation of the supercritical instability for charge impurities in graphene remained elusive until 
recently. One can reach the supercritical regime by collecting a large enou^ number of charged adatoms in a certain 
region of graphene. Such an approach was recently successfully realized [y| by using a STM tip in order to create 
clusters of charged calcium dimers. 

The stu^ of supercritical instability of one Goulomb center in the continuum model in gapped graphene was 
extended [3] to the case of the simplest cluster of two equally charged impurities when the charges of impurities 
are subcritical, whereas their total charge exceeds a critical one. We determined the critical distance between the 
impurities separating the supercritical and subcritical regimes as a function of the charges of impurities and th e g ap. 

An interesting problem of two oppositely charged impurities in graphene was recently considered in Refs.(8|, l9|. 
By studying the Dirac equation for quasiparticles in graphene in the field of point electric dipole, it was shown that 
this equation admits towers of infinitely many bound states exhibiting a universal Efimov-like scaling hierarchy. The 
electric dipole problem is obviously particle-hole symmetric one and, consequently, the energy levels are symmetric 
with respect to the change of the sign of energy E —>■ —E. Therefore, naively the bound electron and hole states can 
cross as the dipole moment increases only at A = 0. However, such a level crossing is impossible for two states with 
the same quantum numbers in view of the avoided crossing theorem 0- Indeed, the explicit calculations in Ref.Q 
showed that the positive and negative energy levels first approach each other and then go away. Since no electron 
bound state dives into the lower continuum as the electric dipole moment increases, the conclusion was made in Ref.Q 
that supercriticality is unlikely to occur in the electric dipole problem. 

In this paper, we address the problem of supercritical instability for quasiparticles in graphene with two oppositely 
charged impurities situated at finite distance R {finite electric dipole). We show that a new type of instability is 
revealed in this problem connected with the change of localization (migration) of the electron wave function on 
impurities as the distance between the two impurities changes. This migration leads to the spontaneous creation of 
an electron-hole pair with the electron and hole partially screening the positively and negatively charged impurities. 
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respectively. The necessary condition for the instability to occur is the crossing of the energy levels of the corresponding 
single positively and negatively charged impurities. 


II. DIRAC EQUATION 


The electron quasiparticle states in the vicinity of the K± points of graphene in the field of two oppositely charged 
impurities are described by the following Dirac Hamiltonian in 2 + 1 dimensions (we set = 1 in what follows): 

H = vf<tp +£,^(Jz - eVir), (2.1) 

where —e < 0 is the electron charge, p is the two-dimensional canonical momentum, ai are the Pauli matrices, and 
A is a quasiparticle gap. The gap in graphene can be opened in various ways, e.g., by depositing graphene on a 
substrate [ij, [I 4 I or simply due to finite-size effects in graphene nanoribbons . Hamiltonian (EU acts on two- 
component spinor which carries the valley (^ = ±) and spin (s = ±) indices. We will use the standard convention: 
'I'+s = {'^At'4’b)k+s^ whereas = ['>Pbt'4’a)k-s, and A,B refer to two sublattices of hexagonal graphene lattice. 
The regularized interaction dipole potential for charged impurities iQ, Q = Ze, situated at (±i?/2,0) in the {x,y) 
plane is given by 


V (r) 


Q f 1 

K (x + RjZf- +y'^ +rl 


{R^-R)] 


( 2 . 2 ) 


where k is the dielectric constant and the Coulomb potential of each impurity is regularized by tq, which is of the 
order of the graphene lattice spacing. Since the interaction potential ( 12 . 21 ) does not depend on spin, we will omit the 
spin index s in wave functions in what follows. Further, for the sake of definiteness, we will consider electrons only 
in the valley (the Dirac equation for the electrons in the A"_ valley is obtained by replacing A with — A). The 
main difficulty in solving the Dirac equation for the electron in the electric dipole potential is that variables in this 
problem are not separable in any known orthogonal coordinate system. Therefore, we will utilize in this paper the 
technique of linear combination of atomic orbitals (LCAO) and variational Galerkin-Kantorovich (GK) method. 

Hamiltonian (12.11) has a particle-hole symmetry expressed by = —H, where the unitary operator Q = axR-x 

with the operator TZx of reflection x -A —x satisfies = 1. It follows then that an eigenstate 4'£;(a:,?/) with energy 
E has a partner '^_E{x,y) = ^^E{x,y) = ax'^E^—x^y) with energy —E, hence, all solutions to the Dirac equation 
come in pairs with ±E. 

The Dirac Hamiltonian with the electric dipole potential commutes with the operator U = a^KlZy where K is the 
complex conjugation, TZy is the operator of reflection y -A —y. Since = 1, wave functions are split into two classes 
U\'^x) = A|tI>A), where A = ±1. Since the operator U is antilinear, the function |'I'_) related to the function |4'+) by 
means of the phase transformation, l^*-) = i|'I'+), is an eigenfunction with A = —1. Therefore, there is no need to 
consider functions Hence the components of wave functions |4'+) = (^, x)^, satisfy the constraint conditions 

(j)*{—y) = <(>(?/), X*(~y) = ~x{y) which are consistent with the Dirac equation. 

It is convenient to work with dimensionless quantities h = ^ and e = 3 and use dimensionless coordinates and 
distances defined in units of i?A = We introduce also dimensionless coupling constant C = Thus, Dirac 

equation reduces to a system of two coupled ordinary differential equations of the first order 

/ -i{dx + idy)(j) -f (w - e - l)x = 0, , . 

\-i{dx-idy)x+{v-eEl)^ = Q, 


where v{r) = —eU(r)/A. 
wave functions with A = 
LGAO method. 


We will find numerical solutions of these equations by using the GK method in the class of 
1. Still it is instructive to analyse the Dirac equation analytically. For this, we utilize the 


III. LCAO TECHNIQUE 

The LCAO method is well known and widely used in molecular physics H- Wave functions in this method are 
chosen as linear combinations of basis functions, where the latter are usually the electron functions centered on the 
corresponding atoms of the molecule. By minimizing the total energy of the swtem, the coefficients of the linear 
combinations are then determined. The LCAO method was recently used in Ref.[^ to solve the two Coulomb centers 
problem in graphene and the corresponding results are similar to those found by matching of the asymptotics tZl ■ 
However, the LCAO method was not applied to the analysis of the electric dipole problem performed in Refs.[l, QT 
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As to the atomic orbitals for problem (I2.3L we take the wave function of the lowest energy bound state in the 
field of positively charged impurity and the wave function of the highest energy bound state for negatively charged 
impurity (these wave functions are related to each other by charge conjugation). 

We begin our analysis with the Dirac equation for the electron in graphene with one positively charged impurity 
hp'l/p = e'Pp with Hamiltonian 


hp = -i{(Txdx + CTydy) + (Jz 


C 

\/r'^ +rl 


(3.1) 


[The Hamiltonian for the electron in the field of negatively charged impurity is obtained from the Hamiltonian hp 
by the change of the sign of the last term in hp.] 

We determine numerically the energy levels by using the shooting method with regular boundary conditions at r = 0 
for the wave functions and requiring that the wave functions decrease at infinity. The energy of the lowest (highest) 


e 



FIG. 1: (Color online) The energy of the electron bound state with j = 1/2 in the regularized Coulomb potential with ±Q 
charges as a function of ( for different values of the regularization parameter: ro = O.OSAa (red solid lines), ro = O.OIRa (blue 
dashed lines), ro = 0 (green dash-dotted lines). 

electron bound state with the total angular momentum j = 1/2 in the regularized Coulomb potential with the charge 
+Q {—Q) is plotted in Fig[T]for different values of the regularization parameter tq as a function of (. The levels 
which descend from the upper continuum correspond to the positive charge +Q while those which are pushed from 
the lower continuum and grow with ^ correspond to the negative charge —Q. These results are in accordance with 
calculations in Ref.[l^ (see Fig.4 there) where slightly different regularization for the one Coulomb center potential 
was used which admitted an analytical solution. They also reproduce qualitatively the behavior seen directly at the 
tight-binding level on a honeycomb lattice [l^ . For nonregularized Coulomb potential with positive charge the lowest 
bound-state energy is always positive, it reaches the value e = 0 for C = 1/2 and becomes purely imaginary for 
C > 1/2 (the fall into the center phenomenon @-@). For regularized Coulomb potential with charge +Q{—Q), the 
lowest (highest) bound-state energy crosses e = 0 and dives into the lower (upper) continuum at certain value of the 
charge. For example, for tq = 0.05i?A, this happens for Ri 1. 

Since the operator Uc = interchanges the hp and Hamiltonians, UchpU^ = —hn, the electron levels in the 
field of negatively charged center described by the Hamiltonian are obtained by the reflection e —?► — e and intersect 
with the levels of the Hamiltonian /ip at e = 0. The corresponding critical value C,c when this happens will play a 
crucial role in the behavior of energy levels in the dipole potential because the behavior of these levels dramatically 
changes depending on whether C < (/c or </ > C,c- For chosen values of the regularization parameter rg in Fig[Tl the 
critical coupling (c = 0.6 (rg = O.OIRa) and (/c = 0.7 (ro = O.OSRa). In general, (c increases with the increase of rg. 

We are ready now to consider the Dirac equation for quasiparticles in graphene with two oppositely charged 
impurities. The corresponding Hamiltonian has the form 

h =-i(axdx + aydy) + - ^ (3.2) 

V^n + rl + rl 

where . We seek the wave function as a linear combination (hybridization), 

l^-) =z;p|4'p)-f UnlT-n), (3.3) 

of the wave functions ifp and which are eigenstates of the Hamiltonians hp and /i„, respectively, with eigenvalues 
±eo, where eg is the energy of the lowest energy electron bound state in the field of one Coulomb center with the charge 
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+Q. Explicitly, the functions |4/p), I'i'n) with the total angular momentum j = 1/2 are given in polar coordinates by 


'I'P = 


firp) 

-ie^^pg{rp) 


'ir 


-ifirn) J 


(3.4) 


where exp[—= (x ± R/2 — iy)/rp^n- The radial functions f(r),g{r) are computed numerically in the regularized 
potential of one positively charged impurity (see Figl2]). 


f(r),g(r) 



FIG. 2: (Color online) The radial functions in the regularized potential of one positively charged impurity with ( = 0.85, 
ro = 0.05J?a: the function /(r) is shown by blue solid line and g{r) by red dashed line. 


It is substantial for our analysis below to use the electron wave functions I'Fp), j'l'n) in the field of single regularized 
Coulomb centers whose energies may cross zero. By making use of the wave function (13.3L we project the Dirac 
equation h|4') = e|'I') on the states |'I'p) and |4t„) and find the following secular equation: 


det 




= 0 , 


(3.5) 


where hij = {i\h\j), {i\j) = Sij, i,j = 'I'p,'I'„, S = ('I'p|'I'„) = (4'„|'I'p), (4'|4') = 1. It is easy to see that the overlap 
integral, 

S = [ dxdy(f{rp)g[rn )-—— + f{rn)g{rp) ^ ^ ^ 

J \ '^n '^p J 

vanishes after changing x —> —x in the first term in the last equality. In order to calculate the coefficients hy, it is 


e 



FIG. 3: (Color online) The energy of the bound state levels for ro = 0.05i?A as a function of a distance in the LCAO method: 
C, = 0.65 (green dash-dotted lines) and C, = 0.8 (red dashed lines), ^ = 0.9 (blue solid lines). 


convenient to represent the Hamiltonian in the form h = hp + C,!\Jr'^ + rQ = hn — Q/ ^Jvp -|- Tq . Then we find 

c 


Vp “ eo + ('I' 


hpn — (^p 


c , 


r2+r2 


l^p) — Co + CC' — —hr, 


(3.6) 


= -CA = h 


np' 


(3.7) 
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We find that the Coulomb integral C equals 


C' = (p| 


\Ai + ' 


oo 

■,\v) = j 


Ardr 


y^ir + RY + rl 


=K 


ArR 


(r + RY + b 


(/^W +5^(0) 


(3.8) 


where K{k) is the complete elliptic integral of the first kind. We note that the Coulomb integral C is positive definite 
and monotonously decreases with increasing R. Further, the resonance integral A equals 


OO OO 


A = (j>\—^^=\n) = 2 f dx f dyf{rp)g{rn) 

-L I 


X-RI2 


Irf, + 




(3.9) 


The Coulomb integral C and resonance integral A can be computed numerically with the functions / and q found for 
an isolated impurity problem. Asymptotically at large R they behave as C ~ 1/i? and A ~ exp(—-\/l — e^R). 
Finally, we obtain the energy levels 


e — [hppY + {hpnY — i\/ (eq + CC^Y + C^A '^, 


(3.10) 


which are obviously symmetric with respect to the replacement e —> —e in view of the charge conjugation symmetry of 
the problem under consideration. We note that the energy levels never cross in agreement with the avoided crossing 
theorem 0 - Since the Coulomb and resonance integrals, C and A, tend to zero as i? —>■ oo, the energy of the system 
for large distances between the impurities tends to e —>■ ±|eo| as expected. The coefficients of the wave function of the 
negative energy level are given by 


(hnpY T {^pp + \/{hppY + {hpnYY 
hpp ~t“ ^ylJYipY~d~~(JlpnY 
\J{hnpY T {hpp + {hppY + {hpnYY 


(3.11) 

(3.12) 



li'I’O.y) 


FIG. 4: (Color online) The square modulus of the negative energy wave function for — 0.85 and various distances between 
the impurities: R = 1.25 (upper left panel), R = 2.0 (upper right panel), and R — 2.25 (lower panel) in the LCAO method. 

We plot the energy levels of the system for tq = 0.05i?A in Figl3]as functions of i? for C = 0.65, ( = 0.8 and ( = 0.9. 
In the first case, we have C < Cc = 0.7 and the bound state levels monotonously converge to each other as R increases 
and never cross. For the couplings C = 0.8 and () = 0.9 which are larger than Cc = 0.7, their behavior is no longer 
monotonous. For small R, the levels converge like in the previous case. However, after the maximal convergence of 
the levels they go away with the subsequent increase of R. This behavior is typical for the avoided crossing [lOj . 

Ea. (l3.10p and the facts that C and A monotonously depend on R and C > 0 imply that the level repulsion can 
take place only for eg < 0 and the energy levels converge most closely for = |eo| from which we can determine the 
corresponding distance Rm- Exactly at this distance Rm we have Vp = Vn = and, consequently, the probability 

to find the electron near the positively and negatively charged impurities is the same. Furthermore, for C^C k, [cqI, the 
difference of the coefficients Vp and squared equals 

2 2 2/lpp ^hpp + i^ppY ffi i^pnY^^ 

{hnpY + {hpp + y/{hppY + (hpnYY 


sign{hpp) = sign(eo + (C). 


(3.13) 
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The following qualitative picture appears for eg < 0. For small i?, |ti„| ^ |r!p| and, therefore, the electron wave 
function of the negative energy level in the LCAO method is localized mainly on the negatively charged impurity. 
Although this result seems to be counterintuitive, it is quite natural. The point is that the energy spectrum of the 
system for i? = 0 is composed of the upper and lower continua. Since the chemical potential in neutral graphene is 
zero, the electron states of the lower continuum are occupied. For small R, the positively charged impurity produces 
electron bound states which descend from the upper continuum. Obviously, there are also charged conjugated states 
localized near the negatively charged impurity, which rise from the lower continuum as R increases. These states are 
occupied for sufficiently small R. Since Irinj = |np| at the point of maximal convergence R = the probability 
to find the electron near the negatively and positively charged impurities is then equal. As the distance between 
impurities R increases further, the difference of the square moduli (13.131) changes sign because the Coulomb integral 
C decreases with increasing R. This means that the electron wave function changes its localization to the positively 
charged impurity. This change of the wave function localization (the relocalization or ’’migration” effect [ 13 ) is 
explicitly shown in Fig|4]for Q = 0.85 when eg = —0.453 < 0. The value of Rm in general depends on rg and on the 
value of the coupling constant C- As the coupling increases, Rm decreases (see Figl5|). It is clear that if eg > 0, then 


Rm 



FIG. 5: The dependence of the distance Rm on the charge ^ at rg = O.OSAa. 

nothing interesting happens. Indeed, since hpp > 0 for eg > 0, we find that |n„| > |t!p| for any distance R between the 
impurities. Therefore, the wave function is always localized on the negatively charged impurity and, consequently, 
the wave function of the highest occupied state does not change its localization. Note that the behavior of the energy 
levels in this case given by green dash-dotted lines in Figl3]as functions of R is monotonous unlike the case where eg 
is negative. 

The LCAO method used in this section has some limitations. For example, it cannot be applied for small distances 
A, and for values of charges when levels in the field of single centers enter continua. Therefore, we consider below the 
GK variational method which is free of these limitations. 


IV. VARIATIONAL METHOD 

In this section, we use the GK variational method in order to find numerically solutions of the system of equations 
(lOl) . The trial functions which belong to the class A = 1 (see the discussion above Eg. (12.31) 1 and satisfy the asymptotic 
at large distance are taken as 


Hx,y) = A{x,y)\ y J 2 k{x)y -b O J 2 k-i{x)y 


N 








X{x,y) = A{x,y) 


N 

k^O 


N 

E^ 

fc=i 


g 2 k{x)y +> 52fe-i(a;)y 


..2k-l 


(4.1) 


where A{x,y) = exp[—Vl — y (|a;| — R/2)^ + y"^ + Xq]. These functions are substituted into Eg. (12.3D and their 
orthogonality to the residual with respect to the variable y is required. Then we obtain the system of ordinary 
differential equations for the functions fk, gk- 
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The same method can be also used for solving the Dirac equation for the electron in the field of one Coulomb center 
with regularized potential. We found that approximate spectra determined by the GK method for the various number 
of terms in the ansatz N, N describe the exact spectrum in the best way when N = N'. 



FIG. 6: (Color online) The bound state energy levels in the electric dipole potential for = 0.85 and ro = 0.05i?A obtained in 
the LCAO (blue dashed lines) and variational (red solid lines) methods. 

The approximate spectrum of the Dirac equation for the electron in the electric dipole potential is plotted in Fig. 
m for ( = 0.85, ro = 0.05i?A in the N = 1,N =1 approximation (red solid line). For comparison, we plot in the 
same figure the energy levels obtained by using the LCAO technique (blue dashed line). Qualitatively both methods 
give similar results. For ( = 0.85, the energy of the lowest electron bound state in the single positive Coulomb center 
problem (see, FiglT|) is negative (for the chosen charge ( = 0.85 the single positive Coulomb center has only one 
negative energy level). Therefore, the wave function changes its localization on Coulomb impurities (see the video 
in the Supplemental Material 0). Our analysis performed by making use of the GK variational method confirms 
the conclusion made in the LCAO technique that the migration of the wave function of the highest energy occupied 
electron bound state takes place when the charges of impurities are such that the energies of the corresponding single 
Coulomb center problems cross. 

The iV = 1, iV =1 approximation in the GK method was also applied to the study of the spectrum of an isolated 
impurity. A good agreement with the exact spectrum in Fig|T]was obtained. The value of the critical charge of the 
impurities (which form the electric dipole) depends on the regularization parameter. For ro = 0.05i?A, the exact value 
of the critical charge is = 0.71, and the approximate values determined by the GK method are ~ 0.675 

and « 0.68 for the iV = 0,iV =0 and N = 1,N =1 approximations, respectively. 


V. SUMMARY 

In the present paper, we studied the Dirac equation for quasiparticles in gapped graphene with two oppositely 
charged impurities situated at a finite distance (finite electric dipole). By using the LCAO and variational GK 
methods, we found the instability for quasiparticles in graphene in the supercritical electric dipole potential. The 
found instability is connected with the change of localization of electron wave function on impurities and corresponds 
to the spontaneous creation of the electron-hole pair with the electron and hole screening the negatively and positively 
charged impurities, respectively. We showed that the necessary condition for the supercriticality to occur is the crossing 
of the energy levels of single positively and negatively charged impurities. Taken together these two levels should 
traverse the energy distance 2A, which is thus the necessary energetic condition for the supercriticality to occur. 

The instability of the supercritical electric dipole can be observed experimentally by placing oppositely charged 
impurities with C > on graphene and then moving with a STM tip one impurity toward the other one and afterwards 
moving the impurities apart again. The supercritical instability takes place if the impurities become partially screened 
that can be determined by measuring the local density of states. 

We are grateful to V.M. Loktev for useful discussions. This work is supported by the Program of Fundamental 
Research of the Physics and Astronomy Division of the NAS of Ukraine. 
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